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RAPID ITERATIVE REANALYSIS FOR AUTOMATED DESIGN 


By Kumar G. Bhatia * 

Langley Research Center 

SUMMARY 

A method for iterative reanalysis in automated structural design is presented for 
a finite-element analysis using the direct stiffness approach. A basic feature of the 
method is that the generalized stiffness and inertia matrices are expressed as functions 
of structural design parameters, and these generalized matrices are expanded in Taylor 
series about the initial design. Only the linear terms are retained in the expansions. 

The method is approximate because it uses static condensation, modal reduction, and the 
linear Taylor series expansions. The exact linear representation of the expansions of the 
generalized matrices is also described and a basis for the present method is established. 

Results of applications of the present method to the recalculation of the natural fre- 
quencies of two simple platelike structural models are presented and compared with 
results obtained by using a commonly applied analysis procedure used as a reference. 

In general, the results are in good agreement. A comparison of the computer times 
required for the use of the present method and the reference method indicated that the 
•present method required substantially less time for reanalysis. Although the results 
presented are for relatively small-order problems, the present method will become 
more efficient relative to the reference method as the problem size increases. An 
extension of the present method to static reanalysis is described, and a basis for uni- 
fying the static and dynamic reanalysis procedures is presented. 

INTRODUCTION 

Currently, there is considerable interest in automating as much of the structural 
design process as possible. A status report of the current situation in this area is pre- 
sented in reference 1. The papers cited in this reference provide a good summary of 
the work to date in automated design. It is pointed out that the resizing methods for 
structures and the techniques for a systematic way of conducting analyses are two 
important aspects of any automated design system. The present study is concerned 
with systematic reanalysis procedures for automated design. 
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The iterative nature of automated design requires that the individual analyses be 
repeated many times. Therefore, the economy, speed, and accuracy of the individual 
analyses are important factors in the success of an automated design process. A given 
computational technique may thus be well suited for a few analyses, such as in the detailed 
validation of a final design, but the technique may not be well suited for use in an auto- 
mated design process. 

Published work describing iterative reanalysis procedures is heavily oriented 
toward approximate static reanalysis, as evidenced by references 2 to 9. In contrast, 
there appears to be a scarcity of dynamic reanalysis procedures, which are obviously 
more complex and expensive than the static reanalysis procedures. An efficient approx- 
imate dynamic reanalysis is a valuable tool independently and in conjunction with a static 
reanalysis procedure. 

The purpose of this paper is to present an approximate finite -element structural 
reanalysis procedure which is well suited for the iterative requirements of the automated 
design process. It is recognized that the initial formulation of the structural problem for 
iterative reanalysis should be made by taking into consideration the special nature of the 
iterative process, and that such a formulation can be different from the formulation for 
a single analysis. It is shown that the method can be applied to static and dynamic reanal- 
ysis, and a unified approach to static and dynamic reanalysis is presented but the develop- 
ment of the present method in the main body of the report emphasizes dynamic reanalysis. 
In particular, the method as developed herein provides a means of determining the natural 
frequencies and mode shapes of a structure as would normally be required during each 
cycle of an automated design process but without having to perform a complete structural- 
reanalysis. The approximate nature of the method is examined in appendix A. In appen- 
dix B, an extension of the method to static reanalysis is presented, and a basis for unifying 
the static and dynamic reanalysis is described. 

The method may be outlined as follows. The complete structural stiffness and iner- 
tia matrices are assembled by the direct stiffness method for the initial design. A set of 
stiffness design parameters and a set of inertia design parameters are defined so that the 
assembled stiffness matrix is a linear function of the stiffness design parameters and the 
assembled inertia matrix is a linear function of the inertia design parameters. The 
design parameters themselves can be nonlinear (and complicated) functions of the struc- 
tural design variables. In the present study, the design variables are distinct from the 
design parameters and are restricted to structural element thicknesses, cross-sectional 
areas, elastic properties, and/or any combination of these. The complete structural 
matrices are then reduced by static condensation and an eigensolution obtained for the 
reduced system. A modal matrix is then constructed with the desired number of eigen- 
vectors, and the generalized stiffness and inertia matrices for the initial design parame- 
ters are determined. The generalized matrices (for an arbitrary set of design variables) 
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are expanded in two-term Taylor series expansions about the initial design by using the 
design parameters as the independent variables. The analytical expressions to evaluate 
the first derivatives appearing in the Taylor series expansions are derived. Once the 
Taylor series expansions are defined, the generalized stiffness and inertia matrices can 
be directly evaluated for any perturbations in the design variables from the initial design. 
An approximate reanalysis can then be performed by using the generalized matrices 
determined in this manner. 

The method is illustrated for dynamic reanalysis and examples of application of this 
method to simple trapezoidal and rectangular plate models are presented and discussed. 
The results, in terms of natural frequencies, are compared with the results obtained by 
performing a complete structural analysis. A special-purpose digital computer program 
was written to obtain these results. This program is designated RITREAD ^Rapid Iterative 
Reanalysis for Automated Design^. Some of the features of this program are discussed 
in appendix C. The formulation of the method allows the derivatives of the generalized 
structural matrices with respect to design variables to be easily computed, as indicated 
in appendix D. 
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SYMBOLS 

vector of structural design variables, (NV x 1) 
force vector 

function of {d} associated with jKjjj 

system stiffness matrix, (n x n) 

generalized system stiffness matrix, (NM x NM) 

constant matrix of stiffness coefficients associated with ith stiffness design 
parameter, (n x n) 

contribution of the ith element to the system stiffness matrix, (n x n) 

constant matrix of stiffness coefficients corresponding to the jth subelement 
of ith finite element, (n x n) 

system inertia matrix, (n x n) 
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generalized system inertia matrix, (NM x NM) 

constant matrix of inertia coefficients associated with ith inertia design 
parameter, (n x n) 

total number of finite elements 

number of natural vibration modes used in modal reduction 
number of stiffness design parameters 
number of inertia design parameters 

number of degrees of freedom eliminated by static condensation 
number of structural design variables 

number of degrees of freedom retained after static condensation 
total number of degrees of freedom, NZ + NR 

number of constant submatrices required to express jK[j as a linear 
function of design parameters 

stiffness design parameter associated with jKjjJ 

vector of stiffness design parameters, (NP x 1) 

vector of inertia design parameters, (N^ x 1) 

modal matrix, (NZ x NM) 

modal matrix, (n x NM) 

eigenvector corresponding to ith natural mode of vibration, (NZ x 1) 
eigenvector corresponding to ith natural mode of vibration, (NM x 1) 
final transformation matrix, (NR x NM) 
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integer set associated with P m , identifying the subelements for which 
P m =P ij> J = Mi, i = 1,NE 

transformation matrix relating the degrees of freedom eliminated to the 
degrees of freedom retained in static condensation, (NR x NZ) 

vector of structural displacements, (n x 1) 

vector of structural displacements corresponding to degrees of freedom 
retained in static condensation, (NZ x 1) 

vector of structural displacements corresponding to degrees of freedom 
eliminated in static condensation, (NR x 1) 

vector of static displacements, (n x 1) 

square of natural frequency, (rad/sec)^ 

vector of generalized static displacements, (NM x 1) 

column matrix 

rectangular matrix 


Subscript following a parenthesis denotes a partial derivative with respect to a 


design parameter, for example, ([k]). = ([M]) = 

represents a matrix partitioned or separated from a larger matrix. Superscript B 
denotes that the superscripted symbol corresponds to the initial design. Superscript 
is used to denote a matrix transpose. All other symbols are locally defined. 


a[Mj 




Double subscripted matrix 


A TYPICAL APPROACH FOR DYNAMIC ANALYSIS 


A typical approach to dynamic analysis for determining the natural vibration char- 
acteristics of complex structures is described in this section and is used in the subse- 
quent sections as a basis for the development of the reanalysis method. 

In representing a complex structure by an analytical model, a large number of finite 
elements may be needed for a physically continuous structure with an infinite number of 
degrees of freedom to be adequately modeled with a finite number of degrees of freedom. 
The system stiffness matrix [k] and the system inertia matrix [m] can be determined 
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for the complete structure by using the direct stiffness method of reference 10. If a large 
number of finite elements are used for the structural representation, [k] and (m] will 
be of a high order. The dynamic analysis is thus often hampered by difficulties in accu- 
rately and efficiently handling eigenvalue problems of large size. The usual practice is 
to reduce the size of the eigenvalue problem by static condensation in the manner sug- 
gested by Guyan (ref. 11) and Irons (ref. 12). 

The system displacement vector {id} and the system stiffness matrix [k] are 
partitioned as 


and 


(u> 

(nxl) 




( 1 ) 


( 2 ) 


The physical coordinates to be retained are denoted by the subscript h and those to be 
eliminated are denoted by the subscript a. One may, for example, either intuitively 
decide to eliminate some degrees of freedom or use the method of reference 13 to deter- 
mine more rationally the degrees of freedom to be eliminated. However, the method of 
selecting the degrees of freedom to be eliminated is an important practical consideration 
but will not be discussed herein. By the methods of references 11 and 12 a transformation 
matrix [t] relates {U a }- to so that 

(V - [ T K u h} o) 

(NR x 1) 

where 

M - -(Vj'^Kah] (4) 

(NR X NZ) 

The reduced stiffness matrix fk^l is defined by 
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( 5 ) 


[ K h] = [ K hh] + MM 

(NZ x NZ) 


The reduced inertia matrix JM^J can be determined as in references 11 and 12, or the 
same expression can be obtained from the equivalence of virtual work as in reference 14 
(pp. 291-292), and is given by 


N 

(NZ X NZ) 



( 6 ) 


where etc., are the inertia matrix partitions defined by 


[m] = 


[ M hh] 

Kh] 


[Mhot] 

[^aa] 


The reduced structural matrices are used to solve the eigenvalue problem 


| K h] - Ai[Mh]]^i} = {0> 


(i = 1,NZ) (7) 


where ^qA, is the eigenvector defining normal-mode displacements for the NZ degree 

rm 

of freedom system and |[t]^A.j defines the normal- mode displacements for the unre- 
duced system with n degrees of freedom. 


The solution of equation (7) completes the process of obtaining the natural vibration 
frequencies and mode shapes of primary interest. A substantial saving is usually realized 
by reducing the order of the eigenvalue problem from n to NZ as described. The 
order of the reduced structural matrices NZ may still, however, be unwieldy if either 
flutter or dynamic response analyses are to be performed. In order to further reduce 
the order of the stiffness and inertia matrices in such instances, a transformation to a 
smaller number of modal coordinates is effected by means of the familiar normal mode 
approach. A modal matrix [q] is constructed from the first NM eigenvectors corre- 
sponding to the lowest NM frequencies of vibrations obtained from the solutions of 
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equation (7) . The generalized structural matrices are then determined from the following 
equations: 

[kg] = [Q] T [K h ][Q] 

(NM X NM) 

= [Q] T [Khh][Q] + [Q] T [ K ha][s] (8a) 

and 

[MG] = [Q] T [M h ][Q] 

(NM X NM) T 

= [Q] T [Mkk][Q] + fQ] X [M hee ][S] + [[Q] T [M ha ][ S ]] + [S] T [ Mclo ][S] (8b) 

where [sj = [t][qJ. The generalized matrices given by equations (8a) and (8b) are diag- 
onal matrices and are usually of an order much smaller than the system matrices, that is, 
NM « n. The eigenvalue problem using the generalized matrices is of order NM, and 
is expressed by 

[[KG] - ^[MG]]^} = {o} (i = 1,NM) (9a) 

where 

{qi} = [Q]^} (9b) 

In the method for dynamic reanalysis developed in the following section, the order of the 
eigenvalue problem to be solved, if required, is NM. 

DEVELOPMENT OF THE PRESENT METHOD 

The typical dynamic analysis procedure described in the previous section is not well 
suited for efficient reanalysis because (1) evaluation or reassembly of [k] and [m] is 
required for each reanalysis, (2) the transformation matrix [tJ has to be recomputed, 
and [KjJ and [MjJ have to be determined for each reanalysis, and (3) reduction of 
[KjJ and MjJ to [KG] and [MG] is required during each reanalysis. The modal 
matrix [Q determined for the initial analysis could be used as an assumed-mode 
matrix during reanalysis and the order of the eigenvalue problem required to be solved 
for natural-vibration analysis reduced from NZ (as in eq. (7)) to NM (as in eq. (9a)). 

If the natural frequencies and mode shapes of the structure are not required, and only 
[kg] and [MG] are required, the eigenvalue problem of equation (9a) need not be solved 
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during reanalysis. Even then the complete reanalysis procedure as described can be 
very expensive. 

The method presented herein can economically determine the generalized stiffness 
and inertia matrices during reanalysis following their generation during the initial analy- 
sis. The present method involves some modifications in the formulation of system stiff- 
ness and inertia matrices during the initial analysis. The complete system stiffness and 
inertia matrices are expressed during the initial analysis as linear functions of the stiff- 
ness and inertia design parameters, respectively. The generalized matrices and their 
derivatives with respect to the design parameters are then determined. The generalized 
stiffness and inertia matrices are expressed by linear Taylor series expansions in the 
design parameters about the initial design. For subsequent analyses, the generalized 
matrices are determined from the Taylor series expansions. This procedure is devel- 
oped subsequently. 


Formulation of System Structural Matrices 

In the direct stiffness method of finite -element analysis, the system stiffness matrix 
for the complete structure can be generated by properly assembling all the elemental 
stiffness matrices. This assembly is indicated by the relation 

NE 

M ■ ^ pi] < 10 > 

(n x n) i=l 

When only a single analysis is required, there is no further need for jjKjj matrices after 
they have been added to the system stiffness matrix [k], and these elemental matrices 
are not normally saved. When some of the design variables are changed in an automated 
design process, a reanalysis is normally performed. As a consequence of the change in 
design variables, many of the elemental stiffness matrices will have to be recomputed 
before a reanalysis. Considerable computational efficiency can be obtained if the recom- 
putation of the affected elemental matrices can be avoided. This can be done if the ele- 
mental stiffness matrices are expressed in the following form: 

Ki * i f u(<p>)[ K ij] 

(n x n) j=l 
n i 

= I p ij[ K ii] <>» 

j=l 
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where Pjj = fjj is defined as a stiffness design parameter and is a uniquely defined 

function of the design variables {d}; [KjjJ where j = l,nj are constant submatrices 
corresponding to the ith element, and n^ is the number of constant submatrices required 
for the ith element to represent JjCjJ as a linear function of the design parameters. Sub- 
stituting equation (11) into equation (10) yields 

NE H 

M -l I P ij[ K ij] < 12 > 

(n x n) i=l j=l 

The system stiffness matrix as given by equation (12) is a linear function of the stiffness 

o 

design parameters P^. 

The stiffness design parameters are not required to be linear functions of the design 
variables {b}. The concept of the stiffness design parameters and design variables is 
illustrated in the following example for a beam element which includes the effect of shear 
deformation (ref. 15). The stiffness matrix for this element is 


where 

l 





— 

R 

-i R 

-R 

! r 


,2 Ely 

V R+ ~r 

§ R 

,2 Ely 

r R --r 



R 

! r 


Symmetric 


,2 Ely 

t r + -/ 







length of the beam element 
flexural stiffness of beam element 


^ The form of the stiffness matrix given in equation (12) is similar to that published 
by Kavlie and Powell (ref. 4), but was independently developed by the author. 
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R = ( 1 + z3 1 

lk z AG 1 2 E I y j 


k -l 


cross-sectional area of beam element 


shear modulus 


k z constant defining fraction of cross-sectional area of beam effective in shear 

The element matrix may be expressed in the form of equation (11) by writing 
2 

M- 1 p ij[ K «] 

j=l 

where 


p il - EI y 


P i2 = R = '— 


lk z AG 12EI y y 



0 

0 

0 

0 


1 

l 

0 

1 

~ l 



0 

0 

Symmetric 

1 

l 


ii 





The design variables for this beam element could be E, G, Iy, A, b, and h or any 
combination of these variables, and P il and are the design parameters which may 

be nonlinear functions of the design variables . 

Although a large number of finite elements may be required for adequate represen- 
tation of complex structures, some of the design parameters may be identical to 

each other. It would be advantageous to distinguish the number of unique design parame- 
ters from the number of elements as suggested in reference 16. Let jd m , where 
m = 1,NP, be a set associated with a stiffness design parameter P m such that all the 
submatrices for which P m = P^j (where i = 1,NE, j = l,n^ are uniquely identified, 
and there is at least one Pjj = P m . In other words, the various constant submatrices 
associated with any P m (where m = 1,NP) can be summed into a single matrix [K m J 
and the system stiffness matrix expressed in the form 

NP 

M = X P m[ K m] < 13 > 

(n x n) m=l 

where 

M - I I N 

(n x n) ie ^m j e ^m 
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and 

and 


indicates that the summation is performed only for those values of 

j which are identified by the set d m - 
The system inertia matrix can also be expressed as 



i 


N $ 

[M] = ^ (14) 

(n x n) i=i 


where 9i, i = 1,N^, are inertia design parameters, and [MjJ, i = 1,N^, are constructed 
in the same manner as jKjJ, i = 1,NP. The design parameters for the stiffness and 
inertia matrices are defined independently of each other and, in general, will not be the 
same. 

Equations (13) and (14) express the system stiffness and inertia matrices as linear 
functions of two different sets of design parameters {p} and {f}, respectively. The 
recalculation of [k] and [M] for any set of design parameters is facilitated if the con- 
stant matrices and [MjJ are saved during initial analysis for subsequent use. 

In addition to equations (13) and (14), the expressions for the partial derivatives of 
[K] and [m] will be subsequently used to develop the present method. These expres- 
sions are given by 



Equations (13) to (16) will be subsequently used in the derivation of the expressions for 
the generalized stiffness and inertia matrices, and their derivatives with respect to the 
design parameters. 

Determination of Generalized Stiffness and Inertia Matrices During Reanalysis 

The system stiffness and inertia matrices can be expressed as functions of the stiff- 
ness and inertia design parameters (eqs. (13) and (14)), respectively, and the correspond- 
ing generalized matrices (eqs. (8a) and (8b)) can be determined by the method described 
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earlier. The generalized matrices can be expressed in terms of the design parameters, 
and equations (8a) and (8b) are rewritten as 


NP 


and 


[ K °] = I p i + [«] T [ K h a , i ]W 

(NM x NM) i=l L 


[MG] = ' j f> [Q] T [M hM ][Qj + [Q] T rM ha>i ][S] 
(NM x NM) i=l L 


+ [[ < 3j T [ M h a ,l]H 

where jK hh jj, jM hh ^j, etc., are defined by 


T + [s] T [M acl>i ][S] 


( 17 ) 


( 18 ) 


M = 

(n x n) 


and 


>,ij 



'^ho!,iJ 

[W 



■W 


N 

(n x n) 


[ M hh,i] 



>,i] 

[ M ah,i 



M ao!,i 


The generalized matrices are now expanded in Taylor series about the initial design, 
and the series expansions are given by 

NP 

[KG] = [KG B ] + £ ([KG]) (p, - Pj B ) + . . . ( 19 ) 

(NM X NM) j=l J 
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and 


N 9 

[mg] = [mg 8 ] + ^ ([mg]) (^j - 9j B ) + ■ ■ . (20) 

(NM X NM) j=l 

where the superscript B refers to the initial design. 

The following assumptions are now made: 

(1) The generalized matrices are linear in the design parameters. This linearity 
allows neglecting the second and higher order terms in the Taylor series expansions. 
Consequently the generalized derivatives [[KG]]j, j = 1,NP, and [[MG])., j = 1,N^, 
are constants. 

(2) The modal transformation matrix [Q] is a constant matrix. Therefore during 
the reanalysis, the modal reduction (that is, from NZ to NM coordinates) is essen- 
tially the same as in the assumed-mode method except that in the present method the 
assumed modes are the normal modes of vibration (obtained from eq. (7)) of the initial 
structure. It is known that by itself the assumed-mode method during reanalysis would 
give good approximations for the lower modes of the vibration (ref. 17). 

(3) The transformation matrix [t] is a constant for the purpose of calculating 
^[MG]^, i = 1,N^. This assumption simplifies the expression for [[mg|. and will be 

further discussed in appendix A. 

In order to derive an expression for [[KG]j. , equation (17) is differentiated with 
respect to the design parameter Pj. The resulting equation is 

([KGj)j = LQ] T [ K hh,j][Q] + [Q] T [ K ha,j][ s ] + [ c *] T [ K ha]([ s ])j P 1 ) 

(NM X NM) 

where [Q] is a constant matrix. Since [s] = [t][q], 

(NR X NM) 

([»]), = ([t|).[Q] (22) 

From equation (4), 

M > -[KaaJ'VKahl 

(NR X NZ) 
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Therefore 


(M)j = -([K aa ]'‘).[K oh ] - [^''([V]). 

^ J 

or, 

([Tj). = [ K „ a ] ([ K aa]).[ K cm] [ K ah] " ([ K ah]). 

= -[ K «a]" 1 [[ K « a ,i]M + [K ah ,j| (23> 

where jj = ^ K aaj).’ etc- ’ since [ K ] is a linear function of the design parameters. 

Substituting equation (23) into equation (22) yields 

([S])j ■ -[Ka a ] _1 [[K„ a j][s] + [*»!, ,)][«]] < 24 > 

Then 

[Q] T [ K h«](ts]) j = [S] T [K £ , ai j][s] + ts] T [K ohJ ][Q] 

Substituting this expression into equation (21) results in 

([KG]). - [Q] T [K hh ,j]lQ] + [Q] T [K to ,j][S] + [Q] T [K to j][s]| T + (25) 

(NM X NM) 

Equation (25) gives the expression for evaluating the partial derivatives of the generalized 
stiffness matrix. Since the matrix ^[kg]^ calculated from this expression is symmet- 
ric, the generalized stiffness matrix calculated from equation (19) will also be symmetric. 

An expression for the partial derivatives of the inertia matrix can be similarly 
derived by differentiating equation (9) and is given by 
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where ([s]).’ is given by equation (24) . If [t] is assumed to be a constant, ^ 
is simplified to 



([MG]). = [Q] T [M hhj “|[Q] + [Q] T [M h0 ,j][s] ♦ [[Qf[M ha; p] T + [sf[M oa j][8] (27) 

Further discussion of the justification for assuming [t] a constant matrix is presented 
in appendix A. 

Equations (19) and (20), in conjunction with equations (25) and (27), provide the 
basis for an efficient procedure for calculating [kg] and [mg] during an iterative 
reanalysis. The generalized matrices corresponding to the initial design and their 
derivatives with respect to the design parameters can be determined and saved during 
the initial design cycle. For any perturbation from the initial design, the new general- 
ized matrices can be easily determined. 


APPLICATION OF THE PRESENT METHOD 


The results of application and the efficiency of the present method are discussed in 
this section. During the initial analysis cycle, the generalized structural matrices and 
their partial derivatives with respect to the design parameters, are calculated by using 
equations (8a), (8b), (25), and (27). In evaluating the partial derivatives of the general- 
ized inertia matrix by equation (27), the transformation matrix [T j is assumed to be a 
constant matrix. The generalized matrices have been derived by using static condensa- 


3 Equation (26) can be simplified if [t] is assumed to be a constant. It is noted 
that the last four terms in equation (26) contain the product (^Mhajj^ao!] or 
[jVlaa]^ K aaJ > anc * intuitively their sum would be (numerically) substantially smaller 
than the sum of the first four terms. This condition led to the assumption of constant 


for the purpose of evaluating 


([MG]),. 
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tion, modal reduction, and the assumption that they can be expressed by linear Taylor 
series expansions. The generalized matrices during reanalysis are, therefore, approx- 
imate, and the eigensolutions obtained with these matrices are also approximate. In 
order to evaluate the results from the new method, the natural frequencies determined 
from the approximate generalized matrices are compared with the natural frequencies 
determined from a "reference method." 

In the reference method, the system structural matrices are formulated as in equa- 
tions (13) and (14) and are assembled during each reanalysis from the j^KjJ and 
matrices corresponding to the various design parameters. The typical procedure 
described earlier is then employed to determine the eigenproblem expressed by equa- 
tion (7). The natural frequencies obtained from the solution of equation (7) are then 
taken as a basis for the comparison of the natural frequencies obtained from the pres- 
ent method. It is noted that the reference method is also approximate but it is widely 
used; therefore it provides a rational basis for the comparison of the new method. The 
computer central processing unit (CPU) times required for reanalysis by the present 
and reference methods are compared to determine the efficiency of the present method. 

Results of Application 

The present dynamic reanalysis method is applied to simple trapezoidal and rec- 
tangular plates with various boundary conditions. The structure is represented in all 
cases as an assemblage of triangular -plate bending elements (ref. 14, pp. 111-115) 
whereas inertia properties are described by a consistent mass representation. Each 
element has a total of nine degrees of freedom - three normal displacement degrees 
of freedom and six rotational degrees of freedom. The rotational degrees of freedom 
were chosen to be eliminated from the assembled system structural matrices by static 
condensation. 

For the first set of example problems, the trapezoidal plate structure shown in 
figure 1 is represented by 35 grid points defining 48 finite elements. The free structure 
has a total of 105 degrees of freedom. Thickness of each element is taken to be a design 
variable. There are thus 48 design variables and the same number of design parameters 
All percent changes in design variables are with reference to their initial values. The 
various cases studied are described below, and the generalized matrices are 10 x 10 in 
each case unless otherwise mentioned. 

Case 1 .- The grid points 1 to 5 are rigidly fixed and represent a cantilevered bound 
ary condition at X = 0. During the first reanalysis, all the design variables were uni- 
formly increased by 300 percent, and during the second reanalysis they were increased 
by 1500 percent from their respective initial values. The results are given in table I. 

The frequencies from the present method are seen to be identical to the frequencies 
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Figure 1.- Swept, tapered plate, finite -element idealization. Circled numbers are grid points; 
plain numbers are elements. All dimensions are in centimeters (inches). 
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from the reference method; thus, it is- demonstrated that the present method does indeed 
give exact natural frequencies for uniform changes in the design parameters as shown in 
appendix A. 

Case 2 .- The boundary conditions for this case are the same as those for case 1. 

The thicknesses for the elements numbered 2, 3, and 7 were fixed at their initial values, 
and were not changed during reanalysis. The rest of the design variables were increased 
by 25 percent during the first reanalysis, and by 56.25 percent during the second reanaly- 
sis. The results are presented in table II. The maximum error for the first reanalysis 
is about 1.78 percent in the ninth mode, and for the second reanalysis it is about 5.2 per- 
cent in the tenth mode. Thus, the results obtained from the present method can be con- 
sidered to be good even when a few design variables remain unchanged while most of 
the design variables are changed uniformly. 

Case 3 .- The boundary conditions are the same as in the previous two cases. Dur- 
ing the first reanalysis, the thickness of element number 23 was set to zero. In the sec- 
ond reanalysis, the thicknesses for two elements, element numbers 18 and 23, were set to 
zero. From the results in table HI, the maximum error occurs in the second mode; for 
the first reanalysis it is 1.905 percent and for the second reanalysis it is 5.749 percent. 
When three elements were removed, the present method gave an error of about 13 per- 
cent in the second mode and about 14 percent in the tenth mode; these results are not 
tabulated. The results show that the present method gives reasonably good results even 
for drastic changes in the structure. This type of application could possibly be used to 
study the effect of structural cutouts. 

Case 4 .- The boundary conditions are again the same as in the previous three cases. 

A constant concentrated mass of 0.45359 kilogram (1 lbm weight) is placed at the grid point 
number 35. The design variables were uniformly increased by 25 percent and 56.25 per- 
cent of their respective initial values. The results are presented in table IV. The fre- 
quencies calculated from the present method are almost identical to those computed from 
the reference method. Since the concentrated mass remains constant during reanalysis, 
the inertia matrix is effectively experiencing a nonuniform change in the inertia design 
parameters. But the uniform change in the stiffness design parameters causes the trans- 
formation matrix [t] to remain a constant, and therefore the frequencies from the pres- 
ent and reference methods can be expected to be nearly identical. 

Case 5 .- The swept, tapered plate structure is restrained only at grid point num- 
ber 3; a 451.939 newton-meter/radian (4000 lb -in. /radian) spring restrains rotation 
about the X-axis, and the normal displacement in the Z-direction and the rotation about 
the Y-axis are specified as zero. This configuration simulates the arrangement of an 
all-movable control surface. The results presented in table V indicate that the frequen- 
cies from the present method are almost the same as those from the reference method. 
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TABLE HI.- COMPARISON OF NATURAL FREQUENCIES FOR THE CANTILEVERED SWEPT, TAPERED PLATE (CASE 3) 
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TABLE V.- COMPARISON OF NATURAL FREQUENCIES FOR THE SPRING-MOUNTED, SWEPT, TAPERED PLATE (CASE 5) 
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In this example, even though the design parameters are uniformly changed, the transfor- 
mation matrix [t] does not remain a constant because the spring stiffness remains con- 
stant during the reanalysis. The results of this example show that the present method can 
be expected to provide good results for the spring- support type of boundary conditions. 

In a variation of case 5, a concentrated mass of 0.45359 kilogram (1 lbm weight) was 
placed at the grid point number 35. The percent errors in frequencies by the present 
method were almost the same as those for the example without the concentrated mass. 

To study the effect of the number of modes, the example cases 4 and 5 were rerun by 
using the first five modes of natural vibration for the generalized matrices instead of the 
first ten modes of natural vibration. The results showed that there was no appreciable 
increase in the percent errors (five modes compared with ten modes). For the example 
without the concentrated mass, the maximum errors for 56.25 percent uniform change in 
the design variables was 1.096 when using five modes compared with 0.960 when using ten 
modes. 

In the second series of examples, a rectangular -plate structure idealized into 
48 triangular finite elements was considered. This structure is illustrated in figure 2. 

Two types of boundary conditions were considered and these are described by the follow- 
ing two cases. In all the examples discussed below, the generalized matrices were cal- 
culated from the first seven modes of natural vibration. 

Case 6 .- The rectangular plate is rigidly supported along its four edges. The thick- 
nesses were uniformly increased for all the elements by 25 percent and then by 56.25 per- 
cent of their respective initial values. The natural frequencies calculated from the pres- 
ent method were identical to those calculated from the reference method, and therefore 
the results for this method are not presented. In a subsequent study, the thicknesses for 
elements 17 to 32 were increased by 25 percent and 56.25 percent from their respective 
initial values whereas the remainder of the element thicknesses were unchanged from 
their respective initial values. The results for this example are presented in table VI. 

The maximum error in frequency is 2.265 percent in the third mode for the 25 percent 
increase, and about 10 percent in the seventh mode for the 56.25 percent increase. The 
results are considered to be acceptable for reanalysis during automated design. 

Case 7 .- The rectangular plate of figure 2 is simply supported along its edges. For 
uniform increase in thickness for all the elements, the natural frequencies calculated 
from the present method were identical to those from the reference method. The results 
for a case when the thicknesses for element numbers 17 to 32 are increased, are given 
in table VII. The percent errors for this case are somewhat smaller than those for the 
plate with the fixed edges. 
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TABLE VII.- COMPARISON OF NATURAL FREQUENCIES FOR THE RECTANGULAR PLATE WITH 
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Comparison of Computer Times for Reanalysis 


The computer cost saved by using the present method for iterative reanalysis is a 
complicated function of the number of (a) degrees of freedom n, (b) design parameters 
NP and N.^, (c) degrees of freedom eliminated NR, (d) generalized coordinates NM, 
and (e) total number of reanalysis cycles. Table Vm presents a comparison of the com- 
puter CPU times required with the present and reference methods. The computer time 
listed for the eigensolution in the reference method is for complete solution of equa- 
tion (7), and all the NZ eigenvalues and eigenvectors were calculated even though only 
the first NM eigenvectors were needed. The computer subroutine used for the eigen- 
solution is based on Jacobi’s method and finds the solution for all NZ modes. However, 
for the cases in table VIII, the reference method CPU times are not heavily penalized 
since most of the methods that would be used for partial modal solution would still 


require determination of the symmetric matrix 
- 1/2 


■l/2 nT 


K]‘ 1/2 [ K h][Mh 


- 1/2 


or 


jk^j |M h jk^J ' . Thus, the CPU times in table Vin can be used for the com- 


parison and show the efficiency of the present method. Each "analysis" (table VH3) for 
the present method involved evaluating the generalized matrices for the new set of design 
variables and solving the eigenvalue problem of order NM. It is seen that the present 
method yields substantial time savings and hence cost savings when a large number of 
analyses are required. For example, from the first entry in table VIE, the present 
method would require about 3.6 percent CPU time for the reference method if 100 reanal- 
yses are required, and about 22 percent CPU time if 10 analyses are required. This 
statement assumes that generalized derivatives are calculated only once. If the gener- 
alized derivatives are required to be recomputed, the additional time required for each 
recomputation is of the order of the time required for one analysis by the reference 
method. 


A comparison of the number of multiplications required in the present method to 
the number required in the reference method indicates that as the problem size and/or 
the number of reanalysis cycles increase, the computer cost savings achieved by using 
the present method relative to the reference method increase. 

CONCLUDING REMARKS 


A method for iterative reanalysis in automated structural design has been presented 
for a finite-element analysis using the direct stiffness approach. A basic feature of the 
method is that the generalized stiffness and inertia matrices are expressed as functions 
of structural design parameters, and these generalized matrices are expanded in Taylor 
series about the initial design. Only the linear terms are considered in the expansions. 
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The method is approximate because it uses static condensation, modal reduction, and the 
linear Taylor series expansions. The exact linear representation of the. expansions of 
the generalized matrices is also described, and a basis for the present method is estab- 
lished. The formulation of the present method provides a simple means of evaluating the 
derivatives of the structural matrices with respect to structural design variables. 

The method has been illustrated for dynamic reanalysis by applications to simple 
trapezoidal and rectangular plate models. The results in terms of natural frequencies 
were compared with the results obtained from a complete commonly applied structural 
analysis. Several support boundary conditions were considered, and the results obtained 
from the present method were generally in good agreement with the results from the ref- 
erence method. A comparison of the computer times required for the use of the present 
method and the reference method indicated that the present method required substantially 
less time for reanalysis. The present method will become more efficient relative to the 
reference method as the problem size increases. 

The application of the present method to static reanalysis has been described and 
the basic equations have been derived. A procedure for combining the static and dynamic 
reanalysis into a unified approach is offered, but this procedure has not been demonstrated 
by actual application. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., July 9, 1973. 
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APPENDIX A 


THE APPROXIMATE AND EXACT LINEAR REPRESENTATION 
OF GENERALIZED MATRICES 

It was assumed that the generalized stiffness and inertia matrices are linear func- 
tions of their respective design parameters, and that the transformation matrix [t] is 
a constant for the purpose of deriving expressions for the derivative of the generalized 
inertia matrix with respect to an inertia design parameter. It is shown in this appendix 
that [tJ remains a constant if all the stiffness design parameters are uniformly changed, 
that is, the design parameters are changed by a constant percentage from their respec- 
tive initial values. The exact linear representation of the generalized matrices is then 
discussed, and it is shown that the assumption of a constant [t] in deriving the expres- 
sion for the generalized inertia matrix derivatives is justified. 

The transformation matrix [t] for the initial analysis is given by equation (3) and 
is 

M = -[K aa ]' 1 [K a ,h] 

- {*!%,]] l [ P i%hJ 

(AD 

where an arbitrary initial stiffness design parameter P^ B is factored out from [Kq,^ 
and and each element of ^or is obtained by dividing the correspond- 
ing element of ^or by Pj®. Since [^K QIQ! ~| and are linear in 

{p}, each element of and j^ K ahj is either a constant multiplied by the ratio of 

a design parameter to P^® or a summation of several such terms. After a uniform 
percent change in all the stiffness design parameters from their initial values, the ratio 
of any design parameter to the ith design parameter will be the same as the correspond- 
ing ratio at the initial design, and [t] will remain unchanged.^ In this case, the expres- 
sion for ^[MG]j. in equation (27) is exact. In addition, [kg] and [mg] are linear in 
design parameters, and can be exactly represented by the constant and the linear term of 
the Taylor series expansion. 

^ This relationship was pointed out to the author by William C. Walton, Jr., and 
Jerrold M. Housner of NASA Langley Research Center. 


33 
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For a nonuniform change in the stiffness design parameters, the transformation 
matrix [t] does not remain a constant, and the generalized matrices obtained from the 
two-term Taylor Series expansions. are not exact. An exact linear representation of the 
generalized matrices can be obtained if the static condensation is not performed, and the 
eigenvalue problem associated with the complete system matrices [k] and [Mj is 
solved to determine the modal transformation matrix |Q n J. Then 

[kg] = [Q n ] T [K] [Q„] (A2) 

(NM x NM) (NM X n) (n X n) (n x NM) 

and 



(NM x NM) (NM X n) (n x n) (n x NM) 

The generalized matrices thus obtained are linear with respect to the design parameters, 
and the Taylor series expansions obtained by using only the first two terms, are exact. 
The exact series representation is, however, achieved at the cost of solving a considera- 
bly larger eigenproblem , specifically one of order n. As an alternative, one may per- 
form the static condensation and solve equation (7) of the order NZ to form the 
(NZ x NM) modal matrix [q]. Since the transformation matrix [t] has already been 
determined, the modal matrix for the unreduced system Q n J can be calculated. The 
matrix jQ n ~j determined by this method would be, strictly speaking, an assumed-mode 
matrix approximating the exact modal matrix. In any case, can be used to deter- 

mine the generalized matrices which are linear in the design parameters. The general- 
ized matrices can be expressed as 
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and 



where [s] = [tJ[q] 


as before, and 


Q„] = 


[Q] 

[s] 


Equations (A4) and (A5) are identical to equations (8a) (for [KG]) and equation (8b) ]for. 
[MG]), respectively. Differentiating equations (A4) and (A5) with respect to the design 
parameters Pj and , respectively, and assuming j^Q n J to be a constant matrix gives 
expressions for ^[KG]^ and These expressions are identical to equations (25) 

and (27). It should be recalled that in deriving equation (27) for ^[MG]j. , [t] was 

assumed to be a constant. Therefore, the assumption of a constant [t] for the gener- 


alized inertia matrix derivatives is equivalent to obtaining an assumed-mode matrix 
from the normal- mode matrix [q] for the initial design. 


Q 


»] 


/ 
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APPENDIX B 


STATIC REANALYSIS 

The application of the present method to static reanalysis is described in this appen- 
dix. The system of equations required to be solved for the structural displacements is 

[k] {A} = {F} (Bl) 

where {a} is the displacement vector, and {f)> is the applied-force vector. The 
stiffness matrix [k] is formulated as a linear function of the stiffness design parame- 
ters Pp i = 1,NP, as in the dynamic reanalysis. Equation (Bl) can be written in par- 
titioned form as 


[ K n] 

1 

1 

1 

1 

1 

1 

1 

K 12] 

[**i] 

1 

y 22 ] 

l : 



r ~\ 

( A i> 

| > 

( A 2> 

V J 


( F l) 

— ) 

( F 2> ! 

L j 


(B2) 


Let {Ai} be the set of displacements which are to be retained, and be the set 

of displacements which are to be eliminated. For example, could be the vector 

of normal displacements for a wing structure whereas ^ 2 } could be the vector of rota- 
tional displacements. In general, the force vector ^ 2 } corresponding to ^ 2 }) will 
not be identically zero. 


From equation (B2), 
follows : 


^ 2 } can t> e expressed in terms of 


and {^ 2 } as 


( A 2) = [ T s]{ A l} + [ K 22] 1 { F 2} 


(B3) 


where 


[ T s] - -[ K 22]' 1 [ K 2l] (B4) 

Equation (Bl) can then be written as 

[k]{Aj} = <F> (B5) 
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where 



(B6) 


and 

{f> = {Fj} + [t s ] T {F 2 } (B7) 

The order of equation (B5) is less than the order of original displacement equa- 
tion (Bl). For the purpose of reanalysis, [k and {f} can be expanded in Taylor 
series about the initial design. If only the first-order terms are retained, that is, if 
[k] and {f} are assumed to be linear in P^, then 


NP 

[*]-[»]♦ I (fe^-P^ 

i=l 

and 

NP 

{f} = {f b > + £ ({F}).(Pi - P( B ) 
i=l 


(B8) 


(B9) 


where the superscript B denotes the initial design point, P^ are the stiffness design 
parameters, NP is the number of stiffness design parameters, and a subscript follow- 
ing a parenthesis denotes a partial derivative with respect to the design parameter corre- 
sponding to the subscript. The partial derivatives are given by the following expressions: 

(M)i = ([ K ll]). + ([ K 12 ]).[ T sJ + [ K 12]([ T s]). < B10 > 



(B11) 


and 




(B12) 
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During reanalysis, [k] and {¥} can be evaluated from equations (B8) and (B9). 
As a consequence of assuming [fc] and {f} to be linear, the partial derivatives ([k]^ 
and i = 1,NP, are constants and need to be calculated only at the initial design. 

The displacements {A^ can be solved from equation (B5), and ^ 2 } can be deter- 
mined from equation (B3). 

In order to obtain .(A 2 } from equation (B3) during reanalysis, matrices [TgJ and 
[ K 22j * are re( l uire d- Several alternatives are described to obtain these matrices dur- 
ing reanalysis. By virtue of equation (B II) and the fact that the {f} are taken to be 
constants, pTgJ is linear with respect to {p} , and is given by the following expression: 


[ T s] - M + j ([ T s]).( p i - p i B ) 


(B13) 


Note that £k 22 J 1 is not linear in {p}, but may be approximated by assuming it to be 
linear with respect to 1/Pj, i = 1,NP, in a manner similar to reference 7. Thus 


r K r'-fe BT'/f W'/i 1 
M - [_ k 22 J * 1 —jprk - 

i= 1 ( p . ) \ i 

\ 1/ 


■ [ K 22 B ]" + I [ K 22 B ]' ([ K 22]) [ k 22 B ]' 

i=l L 1 


-i|Pi B (Pi B - P t ) 


(B14) 


The matrix triple product [^22^] ([ K 22]j. K 22^| is a constant matrix and has to be 

evaluated only during the initial analysis. It is apparent from equation (B3) that if {F 2 } 
is identically zero, jK^J will not be required. If {F 2 } is small in magnitude rela- 
tive to then jK 22 j * could be assumed to be constant for the purpose of deter- 

mining {A 2 }. Alternatively, it may be argued that assuming |TgJ to be linear in the 
stiffness design parameters implies that [jC 22 J a cons tant. This condition is evi- 
dent from equation (B4) and the fact that is formulated as a linear function of the 

stiffness design parameters. Therefore, assuming [K^J” ^ to be a constant may be 
satisfactory in practice. 

The procedure outlined above - static condensation followed by Taylor series 
expansion - is probably more efficient than reassembling a new [k] matrix for each 
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reanalysis and solving equation (Bl) for {a}. The procedure is, however, not as effi- 
cient as the method for dynamic reanalysis presented in the main body of this report 
because the static reanalysis does not incorporate modal reduction. In order to achieve 
the same degree of efficiency for the static reanalysis as for the dynamic reanalysis, it 
is proposed that the static reanalysis be performed by using an approach parallel to that 
of the dynamic reanalysis, in which case the static and dynamic reanalyses can be inte- 
grated into one unified approach. A procedure to accomplish this integration is described. 

The degrees of freedom to be eliminated by static condensation in the static reanal- 
ysis are selected to be the same as in the dynamic reanalysis. The transformation matrix 
[T S ] will now be identical to the transformation matrix [_T] for the dynamic reanalysis, 
and {A^ and {^ 2 } will be identical to {Uh} and {U®, respectively. If the modal 
transformation matrix [Q] determined for the dynamic reanalysis is used to transform 
the physical displacements {Uh} to a set °* generalized displacements {£}, then the 
generalized matrix [KG] is the same for the static and dynamic reanalyses. The gen- 
eralized displacements ® can be obtained from 

[KG]® = {F g } (B15) 

where 


( F g> ■ raVh} + m t <m 


(B16) 


and 


[s] = [t][Q] 


The vector 




is assumed to be expressed by a two-term Taylor series expansion 


NP 

{ F g}=l F g B } + 1 (( F g» i ( p i - p i B ) 


(B17) 


where 

(c f g», - 


and ([S]). is given by equation (24). If 
constant. 


{F® is identically zero, then {Tg} is a 
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Equation (B 15) can be solved to obtain and -(U^} and (Uq,^ can be deter- 

mined by back substitution; thus 

( u h> = [«]{«} (ms) 

and 

{U«} = [S] {5} + [K a „] ■ 1 {F 0 } (B 1 9) 

Then [s] is determined during reanalysis from 

NP 

[S] = [S B ] + £ ([S]).(P. - (B20) 

i=l 


where ([SjJ. , given by equation (24), is determined only during the initial analysis. 

[Kq-o]" 1 can be determined in the same manner as in equation (B 14), and 

is given by 


M' 1 - 


k; 


b 

a a 


-1 


NP f _i 


3 

''a a 


-i - 1 


p i B ( p i B - p i 


(B21) 


The procedure described for the approximate static reanalysis has not been pro- 
gramed and evaluated. However, the rationale for the static case is similar to that for 
the dynamic reanalysis which is shown to be economical and does give fairly accurate 
results. 
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COMPUTER PROGRAM RITREAD 

The special purpose computer program RITREAD (Rapid Iterative Reanalysis for 
Automated Design), which was developed for solving the example problems presented in 
the main text, is described in this appendix. The program was written in Fortran IV for 
use on the Control Data 6000 series digital computers at the Langley Research Center. A 
simplified flow chart and a brief discussion of some of the salient features of RITREAD 
are included in this appendix. The program listing and the usage procedure are not 
described. 

RITREAD has only a triangular plate bending element and the corresponding consis- 
tent mass representation, and each element has nine degrees of freedom associated with 
it. The finite-element idealizations for the plate configurations of the type shown in fig- 
ures 1 and 2, can be automatically generated by RITREAD from a simple input. 

RITREAD has been written for a maximum of 35 grid points and 48 elements. Each 
element is assigned one stiffness design parameter and one inertia design parameter, and 
each element thickness is taken to be a design variable. Therefore, the number of design 
variables, the number of stiffness design parameters, and the number of inertia design 
parameters are identical to the number of elements. The stiffness design parameter for 
each element is the element thickness cubed, and the inertia design parameter is equal to 
the element thickness. Each grid point has three degrees of freedom - one normal dis- 
placement and two rotations. The rotational degrees of freedom are eliminated by Guyan 
reduction, and the reduced stiffness and inertia matrices are used to solve for the normal 
vibration modes. A maximum of ten modes can be used to construct the modal matrix, 
and the maximum size of the generalized matrices is thus 10 by 10. 

Figure 3 shows a simplified flow diagram of RITREAD. The program is arranged 
in multiple overlays with each overlay performing a logically discernible task. Over- 
lay (0,0) monitors the execution sequence and transfers the required data to different over- 
lays. The large amount of data required for element matrices, and the generalized matri- 
ces are written on the two disk files to save the in-core memory required, and to have 
the capability to eliminate recomputation of the element matrices, the base generalized 
matrices, and the generalized derivatives during a future run for the same structure. 

In overlay (1,0), each unrestrained degree of freedom is assigned an identifying 
number determined from the finite-element definitions which are input into the program 
or automatically generated within the program from a simple input. The constant matri- 
ces of stiffness and inertia coefficients for each element are determined. Each constant 
matrix is n x n, but has only 81 nonzero elements since the plate element has nine degrees 
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Figure 3*- Simplified flow diagram of program RITREAD. 
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APPENDIX C — Continued 


of freedom. These matrix coefficients are written on a disk file as two matrices of size 
9x9 along with a 9 x 1 vector {LG} identifying the element degrees of freedom with 
the system degrees of freedom. For the restrained element displacements, the corre- 
sponding element in {LG} is zero. This procedure minimizes the disk space required, 
and tape 1 can be copied on an external tape and saved for future use, if desired. 

In overlay (2,1), the matrices of the stiffness and inertia coefficients and the asso- 
ciated vector {LG} are read (one element per reading operation) from tape 1. Each 
coefficient matrix is multiplied by the appropriate initial design parameter, and the ele- 
ment stiffness and inertia matrices are thus determined for each element. The contri- 
butions from all elements are added to assemble directly the matrix partitions [Kj^, 
(j^hoj ’ [ M hh]> (j^haj» an< * [ M aa]* These are passed by means of a labeled 

common block to overlay (2,2). 

The transformation matrix [t] is determined in overlay (2,2) by solving the 
matrix equation [K aa ][T] = -|Kj ia J T by use of Cholesky's method. The reduced matri- 
ces [KjJ and [MjJ are then determined by the expressions in equations (5) and (6). 

The modal matrix [Q] is constructed from the solution of the eigenproblem of equa- 
tion (7). The eigenproblem is solved by Jacobi's method, and all the NZ eigenvalues 
and eigenvectors are determined. The matrices [tJ and [Q] are passed to over- 
lay (2,3) by means of a labeled common block. 

The generalized matrices [kG b J and [mG b ] corresponding to the initial design, 
and the derivatives of the generalized matrices with respect to the design parameters are 
calculated in overlay (2,3). The matrix [Sj is first determined from the matrix prod- 
uct [T][Q]. The determination of [kG b ] and [MG B J, and the derivative calculations 
are performed in the same pass through the program. The value of [kG b ] is obtained 
by multiplying the first two terms in the expression for ^[KG]j. in equation (25), by the 
corresponding stiffness design parameter, and summing the resulting matrix over all 
the NP design parameters. Similarly, (jVIG B J is obtained by simply evaluating 

^ ([MG]).^i, where ([MG]). is given by equation (27). 
i=l 1 1 

In determining the derivatives of the generalized matrices by equations (25) and 
(27), the triple matrix products of the form [Qj^jK^h j]M- M T [ K hoy] [s], etc., are 
required. An optimum procedure for obtaining these products is desirable since the 
number of such triple products is 4(NP + N^). The general procedure used in RETREAD 
is illustrated for determining [Q] T |K h j 1 ^[Q]. The 9x9 element stiffness matrix and 
the 9X1 vector {LG} for the ith element are read from tape 1. The element stiffness 
matrix is rearranged into matrix partitions [KZ], QKZRj, and [KRj where [KZj 
is the 3x3 partition for the normal displacements (retained displacements), [KR] is 


43 



APPENDIX C - Concluded 


the 6x6 partition for the rotational displacements (eliminated displacements), and [KZR] 
is the 3x6 partition of the normal displacements - rotations coupling. The vector {Lg} 
is also rearranged so that the first three terms identify the normal displacement degrees 
of freedom corresponding to [KZ], and the last six terms identify the rotational degrees 
of freedom corresponding to [KR]. Any zeros in {LG} indicate that during input the 
corresponding displacements associated with the particular grid points have been speci- 
fied to be zero. The number of nonzeros for the normal displacements IZ and for the 
rotational displacements IR in {LG} are counted. If IZ < 3, the subsequent elements 
to any zero element in the first three elements of {LG} are moved up and the zero ele- 
ments eliminated. If IR < 6, the last six elements of {Lg} are also similarly rear- 
ranged. The rows and columns in the three partitioned matrices, corresponding to the 
zeros, if any, in {lg} are eliminated, and the subsequent rows and columns are moved 
to their proper position to correspond to the newly arranged {LG}. The matrix triple 
products can now be efficiently performed as illustrated by the following expression for 
AM, where [a] = [Q] T [K hhji ][Q] : 

IZ IZ 

A(I,J) = ^ Q(LG(Z),l) £ KZ(*,m)Q(LG(m),j) 

l-l m=l 

Therefore, determination of [Qj^K^h^jO] requires (NM 2 )(lZ 2 ) multiplications 

instead of the (NM 2 )(nZ 2 ) multiplications required if the advantage of the zero ele- 
ments is not taken. The other triple products are also carried out in the manner shown 
for the product [Q] T jK hh 

The matrices [kG®], [mG^J, the design variables {d}, and the generalized stiff- 
ness and inertia derivatives are written on tape 2 in overlay (2,3). The generalized matri- 
ces are calculated for any set of design variables [d) in overlay (3,0). The required 
matrices are read from tape 2, and the generalized stiffness and inertia matrices are 
determined from the expressions for the Taylor series expansion. Once [KG] and 
[MG] are determined, an eigenvalue problem of the order NM is solved for NM 
eigenvalues and eigenvectors. Tape 2 can be copied on an external tape if the same 
problem with different design variables is required to be solved at a future date. 

This brief program description outlines the computational procedure used to imple- 
ment the iterative reanalysis method presented in this paper for solution of the sample 
problems. Although this particular program is limited in size and is restricted to only 
one type of finite element, experience with its use has indicated that the basic program 
organization and solution flow is good and should form a basis for the development of a 
more comprehensive program. 
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DERIVATIVES OF STRUCTURAL MATRICES WITH 
RESPECT TO DESIGN VARIABLES 


In automated structural design procedures based on mathematical programing, the 
derivatives of the structural matrices with respect to the design variables are often 
required. It is shown in this appendix that these derivatives can be obtained in a simple 
manner when the structural formulation is based on the method presented in this report. 

The system stiffness and inertia matrices are expressed in equations (13) and (14) 
as linear functions of the design parameter. The design parameters themselves can be 
nonlinear functions of the structural design variables. Therefore the derivatives of the 
system matrices with respect to a design variable are 



(j = 1,NV) (Dl) 


(j = 1,NV) (D2) 


The derivatives of the generalized structural matrices with respect to a design variable 
are obtained by noting that the derivatives of the generalized matrices with respect to 
design parameters are constant. Thus, 


and 



(j = 1,NV) (D3) 


(j = 1,NV) (D4) 


The design parameters are functions of the design variables, and these functional defini- 
tions are known. Therefore, 3P i y / 9Dj, S^i/SDj, etc., are easily determined. The higher 
derivatives can also be computed on the same basis. 
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The equations (Cl) to (C4) are further simplified if the design variable Dj, 
j = 1,NV, is associated with only one or a few of the design parameters {p} and 
In such a case, the summations do not have to be performed over all the design 
parameters. 
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